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Abstract. The primary difficulty encountered when simulating hypersonic flow is the flow normally in- 
cludes strong nonlinear discontinuities. These discontinuities fall into three broad classes: shocks, slip-lines 
and rarefaction waves. Moreover, in the hypersonic flow regime, the chemistry of hot gases plays a vi- 
tal role and can not be neglected. These facts combine to make the numerical treatment of spacecraft 
reentry a most challenging problem. In this work, we develop a class of finite difference schemes that 
accurately resolve discontinuous solutions to spacecraft reentry flow and are simple to incorporate into 
existing spacecraft reentry codes. 


§1. Introduction. The most simple finite difference schemes for the numerical approx- 
imation of partial differential equations are easy to motivate. Basically what one does is 
replace derivatives by finite difference quotients that approximate them. However, this 
recipe does not alway guarantee a stable scheme. To see this, consider the pde: 


(i.i) 


du . du 
— — -f- A “ — 
di dx 


= 0 


u(x,0) = tio(x). 


Let f n denote nAf with n = 0,1,2,... and Xj denote j As with j = ...,—1,0,1,.... 
Let the grid function u” approximate ti(xj,f n ). One possible finite difference scheme to 
approximate the solution to the pde above is: 


9u 


u? + 1 -v? 
R'plac a At ' 


du IA,- 1 1 — TAJ 1 
Replace — with 1 — - — . 

y dx 2Ax 


u?,, - u? 


This gives and explicit time marching scheme 

At 




Applying this scheme to discrete initial data 


we obtain an exact solution 


u] = {g{0 Ax)) n e i0s * 
g(&) = 1 — ■^■Atsin(0). 


g(6) is called the amplification factor. For all ir > |0| > 0 we have that |<j(0)| > 1. There- 
fore, the amplitude of every wave component with ir > |0Ax| > 0 blows up geometrically 
— this scheme is unstable. 
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A desirable feature of the scheme above is the centered difference approximation for 
dfdx. In addition to its second order accuracy, central space differencing schemes re- 
quire only a three point stencil. There are ways to construct linearly stable central space 
differencing finite difference schemes to solve (l.l). For example, 


u? +1 _ 

; 


_ . du 

Replace — with 
at 

du w;:; - u 

Replace — — with — — 

r dx 2Ax 


At 

,n+l 


n+1 

J-l 


defines what is called backward Euler time discretization. This implicit scheme is uncon- 
ditionally stable. Another implicit scheme that is unconditionally stable is the trapezoidal 
rule. It is given by 


du 

Replace — with 
dt 

du 

Replace — with - 
dx 2 


u” +1 - u? 
At ' 


1 ( n ? +1 - n ?_ 1 
2Ax 


+ 


u 


n+l 

j+1 


IX 


n+l' 

J-l 


2Ax 


While these two implidt schemes are stable and therefore convergent for smooth data, 
they both work poorly for problems that have discontinuous solutions. This fact results 
from their large phase errors in the high frequenaes. It is also possible to build an explicit 
conditionally stable centered scheme by the addition of artifiaal viscosity. A second order 
in space and time scheme for (1.1) is the popular Lax-Wendroff scheme 


u 


n+l _ 


n At , 


n 

j+i 


-«i-i) 


1, At 


+ § ( ^ A )>j+i- 2u i + ^ ) * 


This scheme also produces poor results when the solution to (1.1) is not smooth. To see 
this, consider (1.1) with A = 1 and with initial data 


ii(x,0) = ti o (s) = {J 


if |x — 1/2| < 0.10 
otherwise 


(see the left pulse around x = 0.5 in figure 1.1). The exact solution to this problem 
at t = 1.0 is just the pulse translated to the right by 1 unit, (see the right pulse around 
x = 1.5). The small boxes in figure 1.1 depict the numerical results to this problem coming 
from the Lax-Wendroff difference scheme. 
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Figure 1.1. The solid line represents the ini- 
tial data and exact solution at t — 1.0. The 
small boxes indicate the grid values coming 
from the Lax-WendrofF difference scheme. 



Figure 1.2. The solid line represents the ini- 
tial data and exact solution at t — 1.0. The 
small boxes indicate the grid values coming 
from one of our modified fourth order central 
difference schemes. 


The wiggles axe a result of the high frequency phase error inherent to the Lax-Wendroff 
scheme. 

When simulating the flow of air around hypersonic spacecraft on reentry, large gradi- 
ent flows are the rule and not the exception. It should be obvious from the simple example 
above that traditional finite difference schemes can yield less than satisfactory results in 
these flow regimes. A certain amount of care must therefore be given to the design of 
numerical approximation techniques to accurately solve hypersonic flow problems. Elimi- 
nating nonphysical oscillations that are intrinsic to approximations coming from centered 
schemes applied to nonlinear systems of conservation laws is the main focus of this work. 
We modify a class of central space difference schemes to completely eliminate spurious 
numerical oscillations found in the example above. Compare the results coming from one 
of our modified central schemes presented in figure 1.2. Our modified schemes arc intended 
to be simple enough so as to easily lend themselves to reliably approximate solutions to 
the equations governing the reentry of hypersonic spacecraft. 


§2. The modified centered scheme for the scalar law. In this section, we de- 
velop a simple modification to a central differencing scheme applied to nonlinear, scalar, 
conservation laws of the form: 

( 2 . 1 ) = 

An explicit forward Euler in time and central differencing in space finite difference formula 
applied to this problem is written as 

( 2 . 2 ) u ? + ' (/(»?»/»)■ -/(«?-./»)). 
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where Uj is used to indicate the approximation of ti(x,i) at time level t n and spatial 
grid block centered at Xj. u ” +1 / 2 w ill be use d ^ ere t° denote some approximation to 
•u(Xj+x/ 2 ,i n ). For example, if we take 


«?+!/» = 5K+i+“"). 

then the resulting scheme is second order accurate in space and first order in time. One 
could hardly ask for a simpler approximation scheme to solve a partial differential equa- 
tion. Unfortunately, as demonstrated above, this second order in space central differencing 
scheme is strictly unstable. To stabilize the central space differencing approach, we add a 
parameter free nonlinear artificial viscosity term to the right hand side of (2.2). This added 
viscosity term will have the property that in regions where the numerical approximation 
is not at a local extrema, the numerical viscosity term is identically zero. 

Throughout this section, we make the following assumptions concerning the midpoints 
required by the forward Euler difference scheme (2.2): 

(2.3a) For each j, the midpoint Uj +1 / 2 is given by a continuous function of grid values tifc, 
k = j,j±l,j±2,.... 

(2.3b) For each j, the midpoint Uj+ 1/2 is restricted to lie in the interval given by neighboring 
grid values uj and tij+i. 

The minmod function is defined by 


minm n 



inmCMilyDsW*) 

0 


if x • y < 0 
if x • y > 0 ’ 


and we shall use the notation tli" = — u” and ur “ = u” +1 ^ 2 — v™. At each point 

Xj+ 1/2 midway between grid blocks, we consider left and right point values given by the 
formulae 


(2.4) «j +1 /2 = Uj + minmod (ur” , pu L ” ), 

v j+ 1/2 = u j+i + mmmod(uLj +1 ,P^Rj+ 1 ), 

where p > 0 is a finite parameter. 

Remark 1: If (2.4), with any fixed p>l, is applied to a smooth function u(x) that has no 
local extreme points in the interval (*;-i/ 2 >*j+ 3 / 2 )> then for Ax sufficiently small 


( 2 - 5 ) v j + 1/2 — v j+i/2 — “ 7 + 1 / 2 * 

(The parameter p does not depend in any way on the grid size.) 
Our modified forward Euler centered scheme now reads 


( 2 . 6 ) 



22-5 



where /i/(u 1} ii 2 ) is any continuous 2 point monotone flux function that is consistent to 
/(u) in the sense that 


h f (u,u) = /(u); 

see [4]. Referring back to (2.5), note that we should expect generically that 

hf ( v j+i /2 > ^j+i/a) = /( u i+ i/a)> 

so we should expect that generically (2.6) will reduce to the centered scheme given in (2.2). 
The modified scheme can be written as 

(2-7) «? + ‘ = uj - £ (/K +1/3 ) - /(«?_!/,) 

“(Qj+ 1/2 - Qj- 1/2)) > 

where 

Qj+1/2 = /( u j+ 1/2) — k/( u J+l/2» V i+l/2)> 

so as to indicate the particular form of the modification to the centered scheme (2.2). 

The modified forward Euler central difference scheme (2.6) satisfies in an obvious way 
the following theorem: 


THEOREM 1. Suppose all intermediate points tx? +1 ^ 2 satisfy (2.3). Let A t& denote the 
time step stability limit for the first order monotone scheme 


A t 


«" + - = «? - 


Jf the time step size At in (2.6) is taken so that 

At < At fcl 
P+1 

/or the /ixed 1 < p < 00 in (2.4), then we have for each n > 0 

min 1x9 <u^ < maxtx9. 

3 3 

Moreover, if we use Var{uj) to denote the usual pointwise variation of a grid values uj, 
we have for each n > 0 


Var(uj) < Var(tx9). 

Remark 2: From the statement of Theorem 1, one should ask the following paradoxical 
question: Remark 1 seems to imply that if the solution to the pde that (2.6) approximates 
is smooth, then (2.6) reduces generically to the strictly unstable scheme (2.2). How can 
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the claim of the theorem then be correct? This question will be addressed at the end of this 
section. 

Proof of Theorem 1: In each cell Ij = (zy_i/ 2 j a: ;+i/ 2 ) determine a number 6j such 

that 


(1/2 + - «?) + (1/2 - e s ) W j+il2 - «J) = o. 

From the definitions of and t>j +1 y 2 we have for each j 


\0j\ < 1/2 


P-1 
P + 1 * 


Cut Ij into two pieces with the dividing point given by xy = Xj + Oj Ax; that is Ij = 
(xy_i/ 2 ,xy) U [xy, xy +1 / 2 ) = ij U IJ. Redistribute data u n onto a nonuniform mesh given 
by Uy(lj UJJ), and let u n denote the redistributed data: 



if x € ij 
if x€ij* 


Scheme (2.6) amovmts to nothing more than a classic monotone finite difference scheme 
applied to data ti n on a nonuniform mesh, and reaveraged back onto the original mesh. 
Since the minimum grid spacing of the nonuniforxn mesh is 


1 

P+1 


Ax, 


the results in [4] make the results of the present theorem obvious. 


The backward Euler version of (2.6) reads: 


( 2 . 8 ) 



where for the implicit scheme we define 



u n+1 
U S- 1/2 


u 


n+1 

i 



u 


n+l 

i+1/2 



to get 

vj+i /2 = u ] +l + xninmod(uiy +1 ,puEy +1 ), 
Vy+i /2 = ■“”+! + minmo d(tx2 y +i , pURj+i)- 


We assume that the midpoint values Uy^ 2 are given by a continuous function of grid 
values at time level n + 1, and that they satisfy (2.3). Therefore, both Vy +1 / 3 and vJ +1 / 2 
are continuous functions of grid values. 
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THEOREM 2. The results of the previous theorem remain valid for the backward Euler 
scheme (2. 8) without any CFL restrictions regardless of the size of 1 < p < co. 

Proof: The idea here is quite simple. Consider a compact and convex subspace of the 
space of sequences defined by 

V = {x : Var(zj) < Var(ti?)} H {a : min(ti?) < xj < max(u?)}. 

i j j j 

Next, define 

v j+ l/afc) = x i +minmod(£fl J -,pxL i ), 

u J+i/2( s ) = sj+i + minmod(zE i+ 1 ,pxR i+1 ), 

and A Fj(x) by 

AFj(x) = hf{vj+i/2{ a )i v j+i/2 (*)) “ M t, J-i/2(*)> v !-i/2( a! ))- 

We use 

< +1 = u? - •^AF i (u n+1 ). 

J J Ax 

to denote the backward Euler scheme (2.8). We will show: (i) The implicit equation above 
has a solution, (ii) The solution obeys the stability results of Theorem 1 . In order to 
accomplish these goals, consider the continuous map p on V defined by 

Fj(x) = Xj~e (xj + ^AFj(i) - . 

P can be rewritten as 

W*) = (!-«) (*; - + «?• 

Prom Theorem 1 it is clear that for e small enough 

Vox(Fj(x)) < (1 - e)Vax(xj) + eVar(u”), 

(1 — e)min(a:j) + emin(u") < Pj(x) < (1 — c)max(ij-) + 

Therefore, we have that P : V — ► V. So by the Schauder fixed-point theorem [1], there 
exists a xi n+1 € V such that ti n+1 = < 7 r (u n+1 ), and the fixed-point u n+1 is a solution to 
the backward Euler scheme (2.8). Moreover, since ti n+1 € "D it also satisfies the desired 
stability properties. 

A second order central differencing scheme is generated by the intermediate point 
value formula 

(2.9) tij+i/2 = -(ty +1 + Uj). 
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Cleaxly (2.9) satisfies conditions (2.3a) and (2.3b). We now derive a fourth order in space 
centered scheme. Let c, + 1 / 2 (x) denote the cubic polynomial that interpolates the cell av- 
erages Uj-i, Uj, u j+1 and u J+2 on cells (x i _ 3/2 ,x J _ 1/2 ), ( z j-i/ 2 > * 3 + 1 / 2 ), ( *j+i/ 2 ,*j+ 3 / 2 ) 
and ( x j+ 3 / 2 i x j+ 5 / 2 ) respectively. If c ;+1 / 2 ( z ) is applied to the values of a smooth func- 
tion u(x, •), we have that Cj +1 / 2 ( x j+i/ 2 ) = u ( x j+l/ 2 >') + c ( s i+i/ 2 )> where e(x) depends 
smoothly on x and e(x) = 0(Ax 4 ). Therefore, since 

1 r x 3+ i/j d 1 . 

Ax J ^/(«(*,<)) <k = — (/( ti ( I i+ 1 / 2 » < )) “ /(«(*i-l/ 2 ,<))) . 

we have that this is equal to 

^ (/( c j+i/ 2 ( z j+ 1 / 2 )) - /(c,— 1 / 2 ( 23 — 1 / 2 ))) + 0(Ax 4 ). 

The formula for Cj + 1 / 2 (xj +1 / 2 ) is 

( 2 - 10 ) c i+i/2(®j+i/2) = |(gj + 9r), 

where, 

<11 = ^( u j+i + «j) ~ g(«j+i - 2u j +«j-i) 

= ^(uj+i + Uj) “ ^(«j +2 - 2 u i+1 + Uj), 

however, Cj +1 / 2 (xj +1 / 2 ) does not always satisfy (2.3). Verification of (2.3) is crucial! We 
can ensure that (2.3) is satisfied by modifying qi and/or q r when necessary. We have taken 
the following approach: Let 

( 2 . 11 ) q t = ^(uj+i + Uj) 

- | min(|tij + i - 2uj + Uj_i|, 3|u, +1 - tij|)sign(ty +1 - 2u,- + u,-i), 

= ^(Uj+l +tlj) 

- \ min(|u i+2 - 2u i+1 + tij|,3|uj + i - Uj|)sign(u i+2 - 2u i+ i + uj), 

D 

and note that qi and q r verify (2.3). So we define 

(2.12) u i+1/3 = |(gj + q r ) 

and note that (2.12) satisfies (2.3) as well. However, when (2.12) is applied to a smooth 
function away from local extrema, we have that for sufficiently small Ax 

Uj+l/2 = 2 ( 9 * + 9r) = Cj+i/ 2 (Xj+i/2), 


22-9 



therefore giving the desired accuracy away from extreme points of the solution. 

We conclude this section by addressing Remark 2. Consider the linear pde with 1- 
periodic boundary conditions 

(2.i3) 1 “+^“ = ° 

«(0,i) = u(l,i) 
u(s, 0) = sin(27rz). 

To approximate the solution to this pde, we will use the following fourth order in space, 
forward Euler scheme: 


(2.14) 

where 


At 


ti n+1 — ti n — 

u j ~ u j £3. ^‘+1/2 v :-il 2 h 


v j+i/2 = + “inmod(ujRj,pux,”) 

= «y _ 1/2 - u] ur} = u ? +1/2 - u] 

u j+i/2 = + 9r), see (2.11). 

Notice that at grid points where the approximate solution generated by (2.14) is not near 
a local extremum, we have that 


v J+i/2 = + $0, see (2.10), 


and one can easily verify that (2.14), using these point values, is strictly unstable - for every 
nontrivial Fourier mode. However, the results of Theorem 1 imply that the approximate 
solution generated by the full nonlinear scheme (2.14) is stable. What in fact happens 
is the following: The scheme away from extrema is unstable and all Fourier modes are 
geometrically amplified. When an oscillation begins to form, the nonlinear limiting mech- 
anism kicks in where local extreme points appear. This limiting will tend to flatten out 
extreme points producing a staircase effect. This phenomena is demonstrated using (2.14) 
to approximate the solution to (2.13) using Ax = 1/100 and the results are plotted in 
figure 2.1 at t — 1.0 


4 
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Figure 2.1. The solid line represents the exact 
solution to (2.13) at < = 1.0. The circles 
depict the approximate solution coming from 
(2.14). 



Figure 2 J2. The solid line represents the exact 
solution to (2.13) at t = 1.0. The circles 
depict the approximate solution coming from 
(2.15). 


While the results above do not contradict the results from Theorem 1, they should however 
be considered unsatisfactory. To eliminate this staircase phenomena, (this phenomena is 
inherent to our modified forward Euler central schemes), we introduce a Lax-Wendroff 
type of viscosity to the right hand side of (2.14). Specifically, we modify (2.14) to 


(2.15) 

where 


u 


n+l _ 


At 


u ? (( v i+i /2 ~ ^i+ 1 / 2 ) ~ ( v j— 1/2 “ Sj-i/i))i 


Sj+l/ 2 “ 2 [Ax) V j+1 ' 7 ~ Vj -^ 2 ) * 


This modification leads to a scheme that is genetically fourth order in space and second 
order in time. More importantly however, in regions where the solution to the pde is 
smooth, (2.15) reduces to a stable numerical scheme. The results from (2.15) are plotted 
in figure 2.2. Implicit schemes of the form (2.8) do not exhibit the staircase phenomena 
demonstrated above for explicit schemes. The viscosity modification hinted to above in 
(2.15) will be analyzed and generalized in future work. 


§3. Hyperbolic systems and numerical examples. The extension of the implicit 
scheme (2.8) to hyperbolic systems is straight-forward. Numerical results coming from the 
implicit scheme applied to hypersonic reentry problems will appear in a future paper with 
C. P. Li of NASA Johnson Space Center. The extension of the modified explicit scheme 
(2.6) to time dependent hyperbolic systems is not so straight-forward. As demonstrated 
in the e xam ples depicted in figures 2.1 and 2.2, some type of numerical viscosity must be 
appended to the basic central scheme. At present, we take the approach outlined below. 
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Consider a one space dimensional, hyperbolic system of m equations 


(3- 1 ) ^/( u ) = °> 

u € R m , / : R m -♦ R m , 

and let Df(u) denote the Jacobian matrix of /(«). Since (3.1) is hyperbolic, the matrix 
Df{u) has m real eigenvalues, denoted here by Afe(u), k = l,...,m and m independent 
eigenvectors, denoted here by r*(u), k — 1,... ,m. We let R{u) denote the matrix -whose 
columns are the right eigenvectors of Df(u), 

-B( u ) = ( r i(u) ”• rm(u)), 

and 

L(u) = R~ l (u). 

The domain and range of the minmod function is extended to R m component-wise. That 
is, for each 1 < t < m 

(minmod(s,t/))j = minmod(si,t/,-). 

The modification of (2.4) to systems is to perform the limiting in so-called characteristic 
variables. Written in matrix form, this procedure is given by 

(3.2) ^i+i /2 = *j+ I K u 'j) nrinmod(£(ti ? ) uh” > pL(tij )ul; )> 

Vj+ 1/2 = minmod(£(n7 +1 )irLj fl ,pI(u7 +1 )tI^ * +1 ), 

The simplest numerical flux hf(u i,u 2 ) we can envision (see (2.6)) is the Lax-Friedichs 
numerical flux function. It is given by 

(3.3) h$ F (u ~ ((/(uO + /( tt 2 )) - i/(«i “ «2)) > 

where the parameter v is an artificial viscosity satisfying 

v > max(max(jA fc (ui)|, |A*(u 2 )|)). 

B 

Since generically we will have that tii = u 2 , the Lax-Friedrichs viscosity term above will 
vanish so that generically we will have that 

h) F {u i,u 2 ) = /(ui) = /(u 2 ). 

Finally, the viscosity modification used in (2.15) will be extended to nonlinear systems of 
equations in the following way. Define cl; and crj by 

(3.4) cl; = minmod(i(ti" )ul; i )urj), 
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c Rj = mmmod(L(uj)uRj, pL(uJ)uiJ), 

and from these define 

SLj = -{cRj - c£j) {cl j/ CRj) , 

SRj = -(CRj - CLj) ( CRj/cLj ) , 

Note that both sij and s R j are well defined and continuous functions of ul'J and urJ. 
Now in each cell j and for each vector component fc, compute 


(3.5) 


{d L j)k = 1/2— min(A i (ti J -)»°) 2 ( 3 Li)** 

{d Rj )k = 1/2— - max(Ai ! (u J , ),0) 2 (sji J -)i : . 


The full numerical flux function we use in our forward Euler calculations below is given by 


( 3 -«) h n + 1 /2 = W p (sj +1/2 , «; +I/2 ) - (B(»? +1 )a ij+l + RW)d Bj ) , 

giving the scheme 

u l +1 = " As { hf i+ 1 / 2 ~ hf 3~V *) ' 


We consider the fourth order modified central scheme (3.6) taking the limit factor 
p = 3.0 and using a Courant number CFL = 0.25. This scheme is applied to the one 
dimensional Euler equations 


d a 

dt pJr dx 


pu = 


0 , 


8 9 , i „ 

H Xm+-( /ra +P)=0, 

S a , , „ 

a e+ & (e + ?)u = °- 


In the equations above, p is a fluid’s density, u its velocity and e its total energy. The 
pressure p is given by the equation of state p = (7 — l)(e — pu 2 ) where the parameter 
7 = 1.4. The Riemann problem initial datum 


(p(x) } u(x),p(x)) 


f (.445, .698, 3.528) ifx<0.5 
{ (.5,0, .571) if x > 0.5 ’ 


defines what is frequently referred to as the Lax Riemann problem. This datum is inte- 
grated to time t = 0.15. See [7,2,6] for comparisons. 100 equally spaced points on the 
interval (0,1) are used for all results presented here. In all figures below, the solid lines 
represent the exact fluid density and the x's denote its numerical approximation. 
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Figure 3.1 depicts the results of the second order Lax-Wendroff finite difference scheme 
applied to the Lax Riemann test problem above with CFL = 0.8. Figure 3.2 depicts the 
results of the first order Lax-Friedrichs finite difference scheme applied to this same test 
problem. Figure 3.3 depicts our new scheme (3.6) applied to the Lax test problem. 



Figure 3.1. The solid line represents the ex- 
act solution to the Lax Riemann problem. 
The x's are results from the Lax-Wendroff 
scheme. 



Figure 3.2. Results coming from the Lax- 
Friedrichs scheme. 



Figure 3.3. Numerical results coming from 
our fourth order scheme (3.6). 


The next test problem is referred to as the bloat wave problem. This problem tests 
the performance of a numerical scheme on very strong shocks and shock wave interactions. 


22-14 



Its initial data contains two shock waves with pressure ratios of 100,000 and 10,000. The 
initial data to this problem is 


(p(x),u(a),p(x)) 


'(1.0,0.0,1000.0) if i <0.1 
<(1.0,0.0,0.01) if 0.1 < z < 0.9 , 
(1.0,0.0,100.0) if z > 0.9 


The boundary conditions at z = 0 and z = 1 are reflecting. The Lax-Wendroff scheme 
will not solve this problem. The initial shocks are so strong that negative densities and 
pressures occur almost immediately. The Lax-Friedrichs scheme will solve it - be it not 
very well. Figure 3.4 depicts the Lax-Friedrichs results using 400 points comparing it to 
the “exact” solution which came from (3.6) using 1200 points. (The exact solution to this 
problem can not be found in closed form.) Figure 3.5 demonstrates the performance of our 
scheme (3.6) on the blast wave test problem again using 400 points. We should point out 
that no modifications of the scheme was required to avoid negative densities and pressures. 




plied to the blast wave test problem. 
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